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Abstract 

We study the effects of noise-intensity fluctuations on the stationary and dynamical properties 
of an overdamped Langevin model with a bistable potential and external periodical driving force. 
We calculated the stationary distributions, mean-first passage time (MFPT) and the spectral am- 
plification factor using a complete set expansion (CSE) technique. We found resonant activation 
(RA) and stochastic resonance (SR) phenomena in the system under investigation. Moreover, the 
strength of RA and SR phenomena exhibit non-monotonic behavior and their trade-off relation as 
a function of the squared variation coefficient of the noise-intensity process. The reliability of CSE 
is verified with Monte Carlo simulations. 
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I. INTRODUCTION 



Langevin models have become increasingly important in modeling systems subject to 
fluctuations. These models have a wide range of applications in physics, chemistry, electron- 
ics, biology, and financial market analysis. In many applications, fluctuations are modeled 
in terms of white noise, which has a delta function correlation with constant noise intensity. 
In general, fluctuations are space-time dependent phenomena; hence, the noise intensity 
fluctuates temporally and/or spatially. Nevertheless, white noise has been be used to model 
fluctuations because at the typical level of physical description, variations in noise intensity 
can be ignored. However, if the variation in the noise intensity fluctuations is large and 
if it occurs in time scales comparable to the physical description of interest, the effects of 
such fluctuations have to be taken into account. Noise intensity fluctuations due to en- 
vironmental variations are particularly important in biological applications. For instance, 
the stochasticity of a gene expression mechanism is derived from intrinsic (discreteness of 
particle number) and extrinsic (noise sources external to the system) fluctuations. Because 
extrinsic fluctuations are subject to biological rhythms with different time scales []], their 
noise intensity varies temporally. 

In financial market analysis, stochastic volatility models (e.g., the Hull & White model 



and Heston model) incorporate temporal noise intensity fluctuations 



51. The stochastic 



volatility mode 
superstatistics 



s assume that noise variance is governed by stochastic processes. In physics, 



ll| take spatial and/or temporal environmental fluctuations into account. 



Superstatistics has been applied to stochastic processes, and it has introduced noise inten- 



sity fluctuations 
A previous study 
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4l6|. by calculating stationary distributions in a Bayesian manner. 
17] indicates the similarity between distributions of a stochastic volatil- 
ity model and Tsallis statistics, which has the same stationary distribution (g-Gaussian 
distribution) as superstatistics in specific cases. 

Most discussions on stochastic volatility models are limited to linear drift terms; hence, 
the application of such models to physical, chemical, or biological systems, accompanied 
by nontrivial drift terms and multiplicative noise, is nontrivial. In our previous paper 18], 
we proposed an approximation scheme that can be applied to general drift terms. We 
considered Langevin equations where the white Gaussian noise intensity is governed by the 
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Ornstein-Uhlenbeck process: 

^ = f(x) + sUt), (1) 
ds 

— = _ 7 ( s _ a ) + ^(t), (2) 

where f(x) is a drift term (f(x) = —d x U(x), where U(x) is a potential), 7 is the relax- 
ation rate, and £ x (t), £ g (i) denote white Gaussian noise with the correlation [{£x{t)Cc(t')) = 
2D x 5(t — t') and (^ s {t)^ s {t')) = 2D s 5(t — if)]. In the present paper, we call the term s£ x (t) 
the stochastic intensity noise (SIN) because the noise intensity is governed by a stochastic 



process. In Ref. [18|], we obtained a time evolution equation using adiabatic elimination 
with an eigenfunction expansion [3]. Although the previously developed method [18!] can 
be applied to nonlinear drift terms, its application is limited to 7 ^> 1 [7 is the relaxation 
rate in Eq. (|2J)]. At the same time, we showed that the time evolution equation of P(x; t) is 
a higher order Fokker-Planck equation (FPE) having derivatives of orders higher than two 
rJ . Analytic calculations of dynamical quantities such as mean- first passage time (MFPT) 
and stochastic resonance (SR) are mainly developed for one- variable FPEs; hence, their use 
in higher-order FPEs is nontrivial. Accordingly, in this paper, we investigate the dynamical 
properties of the coupled equations (OQ) and (j2J), expanding functions of interest (stationary 
distributions and eigenfunctions) in terms of an orthonormal complete set. This technique 
is extensively used to solve FPEs numerically (e.g., the matrix continued fraction method. 
For details, please see Ref. I20] and the references therein). Complete set expansion (CSE) 
can be applied to polynomial drift terms, and it can, in principle, solve for the entire range 



of 7; on the other hand, the adiabatic elimination based method is limited to 7 3> 1 [18 ]. 

In the present paper, we investigate MFPT and SR with a bistable potential [see Eq. 
As stated above, SIN is particularly important in biological mechanisms. In a zeroth-order 
approximation, many important biological mechanisms, such as neuron and gene expression, 
can be modeled with a bistable potential. MFPT and SR have also been extensively investi- 
gated in such biological mechanisms. In the calculation of MFPT, we show that MFPT, as a 
unction of 7, has a minimum around 7 — 1, which is equivalent to resonant activation (RA) 
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26]. Furthermore, by changing p (the squared variation coefficient of noise intensity 
fluctuations [see Eq. (llip ]). MFPT also has a minimum around p ~ 1. In the calculation of 
SR, we show that the SR effect is smaller for smaller 7, which indicates that the SR effect 
is maximized under white noise. In addition, the spectral amplification factor, as a function 



of p, has a minimum around p ~ 1. These results show that the strength of RA and SR 
effects cannot be maximized simultaneously. All the calculations are performed using CSE, 
whose reliability is evaluated via Monte Carlo (MC) simulations. 

The remainder of this paper is organized as follows. In Sec. [Ill we describe the model 
adopted in this study. In Sec. IHH stationary distributions are calculated using CSE. In 
Sec. IIV| we calculate MFPT, which is approximated by the smallest non-vanishing eigen- 
value. In Sec. [V] we investigate the spectral amplification factor of SR by using the linear 
response approximation. In Sec. IVTl we discuss the effects of noise intensity fluctuations on 
RA and SR. Finally, in Sec. IVIH we conclude the paper. 



II. THE MODEL 

We consider the Langevin equations given by Eqs. ([1]) and (T5]) with the bistable potential 

4 2 

U(X) = X ~ - y, (3) 

i.e., f(x) = x — x 3 . In this paper, we investigate the 7 > case for Eq. (J2j) because s(t) is 
constant (s(t) = s(0)) for 7 = 0, and the resulting SIN is equivalent to conventional white 
Gaussian noise. 

By interpreting Eqs. (CQ) and (T5]) in the Stratonovich sense, a probability density function 
P(x, s; t) of (2, s) at time t is governed by the FPE: 

^P(x,s;t)=Jf P(x,s;t), (4) 

where J??o is an FPE operator defined as 

Sf = Sf x + 7 Sf a , (5) 

with 

se. = -± m + tD.!L, ( 6) 

*.-£(.- a) + (7) 

For the asymptotic case 7 — > 00, we used the adiabatic elimination technique to obtain the 
FPE operator isj 

^o = -J^(*) + q|^ (for 7^00), (8) 




FIG. 1: (Color online) Trajectories of SIN for four parameter values of p (squared variation co- 
efficient): (a) p = 0.01, (b) p = 0.1, (c) p = 1, and (d) p = 100. We varied p while keeping the 
effective intensity Q constant. 

where Q is the effective noise intensity given by 

Q = D x (D s + a 2 ). (9) 

Equation (Q is in agreement with the noise intensity of the correlation function, i.e., 
(s(t)£ x (t)s(t'%(t')} = 2Q5(t - t') (see the Appendix). 

From Eq. (171), the stationary distribution P s t{s) of the intensity- modulating term s is 
given by 

p - M = 7m e * p {-^k {s - a)2 }- (10) 

Here, we introduce the squared variation coefficient of the noise intensity fluctuation for 
later use. The squared variation coefficient p is defined as 

D s , s 

or 

where p denotes the squared ratio between the standard deviation and mean of Eq. (fTOI) . 
similar to the Fano factor. Figure [1] shows some trajectories of SIN with (a) p = 0.01, (b) 
p = 0.1, (c) p = 1, and (d) p = 100. These trajectories have the same effective noise intensity 
Q. As p — > 0, SIN reduces to white Gaussian noise with noise intensity Q = D x a 2 . 

In the present paper, the FPE of Eq. (j4j) is solved using CSE and MC. MC is performed 



by adopting the Eu 



method, see Ref. 



er forward method with time resolution At = 10 4 (for details of the 



III. STATIONARY DISTRIBUTIONS 



We calculate the stationary distributions of the coupled Langevin equations ([T]) and (j2J), 
which have been discussed previously [18| for 7 3> 1. The method adopted in this paper 
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is different from the previous one [18j in the range of 7 (previously [18|], it was limited to 
7 ^> 1). In the following, we first investigate the effects of noise intensity fluctuations on 
the stationary distributions. Then, the calculations of the stationary distributions are used 
for the spectral amplification factor in SR (Sec. IVT) . 

The stationary distribution Pq(x, s) of (x, s) has to satisfy the differential equation: 

^ Po(x,s) =0, (12) 

where Jzfo is an FPE operator defined in Eq. (111) . In order to solve Eq. (1121) . we employ 
CSE, which expands Po{x, s) in terms of an orthonormal complete set. This technique is 



extensively used in stochastic processes (e.g., the matrix continued fraction method 20] ). 
CSE can handle systems with polynomial drift terms and it can, in principle, handle the 
entire range of 7. However, in practical calculations, we are restricted to 7 > 0.3 because of 
numerical instability. Considering the symmetry x — > —x in Jzfo and the relation y?2fc(— x) = 
y2fc(^), the stationary distribution P (x,s) admits the even parity expansion: 

K L 

P (x, s) = tp (x)?po(s) 2J ^ C kt np 2 k(x)4>e(s), (13) 

k=0 £=0 

with 



1/4 



* w =(^k) Vw* ( " )exp (4" 2 )- < 15) 

Here, are expansion coefficients, 77 = ^/l/(2D s )(s — a), H n {z) is the nth Hermite 
polynomial, and ( is a (positive) scaling parameter that affects the convergence of CSE. K 
and L are truncation numbers which provide the precision of the obtained solutions. The 
orthonormality and complete relations read 



dx ip k >(x)ip k (x) = 5k> tk , j dsi/j e (s)ip e (s) = 8ff& (16) 

where 5k>,k is Kronecker's delta function. The term ipo(s)ijje(s) forms eigenfunctions of S£ s 
[Eq. ©], i.e., 

^ s [Ms)M*)) = -^o(s)Ms). (17) 
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After multiplying <p2k'(x)ijjei(s)/ ( l Po(x)4'o(s)) by Eq. (Tl2|) and integrating with respect to x 
and s, we obtain the following linear algebraic equation: 



6k 2 

= C kl [2k- — 



c 2 

+C k - 1>e y/2k(2k - 1 



1 - A(2A; _ 1) + 2C 2 ^ {a 2 + 2D S (e + ~ 



-C fc+ i,^ V(2fc + l)(2fc + 2) - C fc _2,^ v/2fc(2A; - l)(2fc - 2)(2fc - 3) 



+2C k . 1/+2 ( 2 D x D sV /2k(2k -!)(£ + 2){£ + 1) + 2C fe _ M _ 2 C 2 AA ^(2*; - 1)^ - 1) 



+AC k - 1>e - 1 C 2 otD x y/2D s k(2k - 1)£ + 4C k ^ 1/+1 C 2 aD xV /2D s k(2k - !)(£+!). (18) 

Because all coefficients vanish for = (0,0), Co,o can be determined by a normalization 
condition [J ds j dxPo(x,s) = 1] as Co,o = 1- The two-dimensional coefficients Ck,e can be 
cast in the form of one-dimensional coefficients C m by the following one-to-one mapping j^ : 

m = 1 + (1 + L)k + £. (19) 

By using Eq. f fl9|) . C k) e can be transformed into C m with 1 < m < M, where M = (1 + 
i^)(l + L). Eq. (fl8|) can be solved using general linear algebraic solvers. CSE transforms 
the differential equations into linear algebraic equations, which are easier to solve. From 
Eq. (jHJ), the stationary distribution Po(x) of x in the asymptotic case 7 — > 00 is given by 

P (x) = J dsP (x,s) = |exp (for 7 ^ 00), (20) 

where Z is a normalizing constant. 

In calculating stationary distributions using CSE, we have to determine K, L, and (. We 
increase K and L until the stationary distributions converge. Although larger values of K 
and L allow better approximation, we find that using excessively large values numerically 
gives rise to divergent distributions. Fig. |2] shows stationary distributions with different 
parameters: D x = 1, D s = 0.1, and a = 0.1 (Fig. GD(a)); and D x = 1, D s = 1, and a = 0.5 
(Fig. [2(b)). Figs. |2ja) and (b) show stationary distributions calculated using CSE for four 
7 values: 7 = 0.3 (solid line), 7 = 1 (dotted line), 7 = 10 (dot-dashed line), and 7 — > 00 
(dot-dot-dashed line). Although the CSE method is valid, in principle, for the entire range 
of 7, it appears that small values of 7 give rise to numerical instability. Consequently, the 
smallest value used in this paper is 7 = 0.3. For 7 — > 00, we used the asymptotic expression 
given by Eq. (120]) . The stationary distributions of MC simulations were computed for four 



FIG. 2: (Color online) Stationary distributions for systems driven by additive SIN. The lines 
and symbols represent distributions calculated using CSE and MC methods, respectively. The 
parameters are (a) D x = 1, D s = 0.1, a = 0.1 and (b) D x = 1, D s = 1, a = 0.5, with 7 = 0.3 
(solid lines and circles), 1 (dotted lines and squares), and 10 (dot-dashed lines and triangles). 
Furthermore, the dot-dot-dashed lines are given by Eq. (|20p . which corresponds to 7 — > 00, and 
the crosses denote MC results with 7 = 100. For viewability, MC data are plotted only for x < 0. 

7 values: 7 = 0.3 (circles), 7 = 1 (squares), 7 = 10 (triangles), and 7 = 100 (crosses). 
Total 10 6 samples each were calculated for the empirical probability densities. Higher peaks 
emerge at metastable sites for smaller 7. The CSE stationary distribution of 7 — > 00 and 
the MC stationary distribution of 7 = 100 are very close, which supports the result that 
a system driven by SIN reduces to one driven by white Gaussian noise with effective noise 
intensity Q. 



IV. MEAN FIRST PASSAGE TIME 



In order to study the dynamical properties of systems driven by SIN, we calculate MFPT. 
With regard to the stochastic volatility model, an escape problem was investigated for the 
extended Heston volatility model in a cubic potential using MC simulations [28j. Non- 
monotonic phenomena such as noise-enhanced stability (NES) 29l-l32| were reported for this 
model. Another study [33| considered a Langevin system, where the temperature (i.e., noise 
intensity) takes two values in a random dichotomatic manner, indicating the occurrence of 
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an RA |2lM26l| phenomenon. 

First, we investigate two basins of attractors and a separatix that separates them in (x, s) 
space. Without fluctuations, the deterministic dynamics of Eqs. flTJ and ([2]) are given by 

! = /(,), | = -,(.-„). (21) 

Considering the quartic bistable potential f(x) = x — x 3 , Eq. (}2~Tj) has three fixed points: 
(±l,a) (stable points) and (0, a) (a saddle point). Deterministic trajectories of Eq. ( 12T1) are 
given by [34| 

rfs = 7 (a - a) 
dx x — x 3 

Specific trajectories, as a function of x, are obtained by solving Eq. (j22|) : 

s(x) = a + iy|x|- 7 |a; 2 -l| 7/2 , (23) 

where W is an integral constant. Figure [3] shows vector field plots of Eq. (I2TT) for three 
7 cases: (a) 7 = 0.1, (b) 7 = 1, and (c) 7 = 10. In Fig. [3], the dotted line represents 
the separatix. We see that the separatix is x = regardless of 7, which is not the case 
for colored-noise-driven systems (the separatix depends on the time-correlation of colored 
noise) . 



Let t s be MFPT to the separatix (x = 0Y 
can be well approximated by an eigenvalue 



or sufficiently low noise intensity, MFPT r s 



35H38|: 



T -=*=h (24) 

where r is the escape rate and Ai is the smallest non-vanishing eigenvalue of the FPE 
operator Jzf [Eq. (j5J]. Equation (|24|) gives a reliable approximation when the noise intensity 
is sufficiently small and Ai is well separated from the remaining eigenvalues [A n (n > 2)]. 
The eigenvalue problem is represented by the equation 

^ 0( j)(x,s) = -X(j){x,s), (25) 

where A and (j)(x, s) are eigenvalues and eigenf unctions, respectively. To calculate the eigen- 
values, we employ CSE as in the case of stationary distributions. According to the sym- 
metry x — > —x in Jdo; the eigenfunctions <f)(x, s) have even [<f) c (— x,s) = <ft e (x, s)} or odd 
[4>°{— x, s) = — 4>°(x, s)] parity symmetry. The even case expansion is identical to Eq. (fl3|) . 
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(a ) a = 1,7 = 0.1 (b) a = 1, 7 = 1 




(c)a = 1,7= 10 
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FIG. 3: (Color online) Vector field plots of Eq. (f2Tj) with a = 1 and (a) 7 = 0.1, (b) 7 = 1, and 
(c) 7 = 10. There are three fixed points at (±l,a) (stables points) and (0, a) (a saddle point), 
which are denoted by filled circles and crosses, respectively. The dotted line is a separatix, which 
separates the two basins. 

and the odd case admits the following expansion: 

K L 

<f)°{x,s) = Po(x)i)o(s)^2^2C k/ Lp2k+i(x)4>e(s). (26) 

fc=0 1=0 

In the same procedure as that for stationary distributions, the even and odd cases of Eq. fl25|) 
can be reduced to linear algebraic equations. By using CSE, Eq. (|25|) for the odd case is 
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calculated as 



C u U,k + \- ^- 2 (2k + l) 2 - 7 f| - C fc+M ^^v/(2fc + 2)(2fc + 3) 

+c fc _i,V2fc(2fc + 1; 



l-^ + 2C 2 J Dja 2 + 2A ' ' 



^2 * l V 2 

-C^-^v^ + l)2fc(2fc - \){2k - 2) + 2C k - lii+2 ( 2 D X D s y/2k(2k + !)(£+ l)(£ + 2) 
+2C k ^ 2 C, 2 D x D s ^2k{2k + - 1) + 4C fc _ 1 ,,_ 1 C 2 J D..av / 2^A;(2A ; + 1)£ 



+4C fc _ 1 ,, +1 C 2 ^a v / 2 J D s A;(2A; + + 1). 

= -AC M . (27) 

Equation (1251) is now transformed into a linear algebraic eigenvalue problem, which can be 
solved with general linear algebraic eigenvalue solvers. 

In practical calculation of Eq. ( 1271) . we increase K and L until the eigenvalues converge. 
In addition, we carry out MC simulations to verify the reliability of the eigenvalue-based 
approximation. MFPT of MC is calculated from the average of the first passage time (FPT) 
of 20000 escape events. For sufficiently small noise intensity, r s can be approximated by 
MFPT r from —1 to because the MFPT dependence on starting points exists only in a 
narrow boundary layer around the separatix [36]. In MC calculation, the initial value is 
x = —1, and s has a Gaussian distribution Af(a, D s ) with mean a and variance D s . Fig. H] 
shows the MFPT (r s ) dependence on 7 and p; the theoretical results obtained using CSE 
are denoted by lines, and the MC results are denoted by symbols. 

Our model includes four parameters: 7, a, D x , and D s . In our model calculations, we use 
7, p, Q, and D x as the given parameters, where Q and p are defined by Eqs. ([9]) and ( fTTI) . 
respectively. When these four parameters are given, a and D s are uniquely determined as 
a = y/Q/{D x (l + p)} and D s = pQ/{D x (l + p)}. First, we investigate the 7 dependence 
of MFPT with D x = 1, Q = 0.08, and various p values. Fig. 0Ja) shows MFPT as a 
function of 7 with four p values: p = 0.01 (solid line and circles), p = 0.1 (dotted line 
and squares), p = 1 (dot-dashed line and triangles), and p = 100 (dot-dot-dashed line and 
crosses). From Fig. IHa), t s is U-shaped and has a minimum around 7 ~ 1, which can 
be accounted for by an RA effect. The conventional RA phenomenon occurs in a bistable 
potential subject to white noise, where the potential fluctuates owing to time-correlated 
stochastic processes. On the other hand, the RA observed in Fig. @](a) is induced by the 
noise intensity fluctuation. Because MFPT increases with increasing potential wall height 
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or decreasing noise intensity (or vice versa), the effect of noise intensity fluctuation on 
MFPT is similar to that of potential fluctuation. This correspondence can qualitatively 
explain the occurrence of the RA phenomenon in the present model. RA induced by a 
noise intensity fluctuation has been reported previously 3^|; it was realized by the random 
telegraph process. As expected, the p = 0.01 case shows a very small RA effect because the 
noise intensity fluctuation is very weak in this case. For larger p, the RA effect is larger 
because the noise intensity fluctuation increases with p [Fig. [TJ. In contrast, the RA effects 
of p = 1 and p = 100 are nearly similar. Remarkably, the effect of RA for p = 100 is not 
larger than that for p — 1, even though the noise intensity fluctuation is stronger for p = 100 

(Fig.QJc) and ( d ))- 

Next, we calculate the p dependence of MFPT by varying p while keeping the effective 
intensity Q constant. Fig. H^b) shows MFPT as a function of p with four 7 values: 7 = 0.3 
(solid line and circles), 7 = 1 (dotted line and squares), 7 = 10 (dot-dashed line and 
triangles), and 7 = 100 (dot-dot-dashed line and crosses). For 7 = 0.3, r s decreases as a 
function of p. On the other hand, t s has a minimum around p ~ 1 for 7 = 1, 10, and 100 (the 
depth of the minimum is smaller for larger 7). As explained, the RA phenomenon is referred 
to as the existence of the minimum as a function of the relaxation rate. The strength of 
the RA effect can be measured by the magnitude of the minima. In all p cases in Fig. Ufa), 
MFPT is minimum around 7 ~ 1. Therefore, MFPT in Fig. 0(b) with 7 = 1 (dotted line) 
can be identified as the strength of the RA effect as a function of p. This indicates that the 
strength of the RA effect increases with p, up to p ~ 1. A further increase in p does not 
increase the strength of the RA effect. 

In Fig. HI the theoretical results obtained using CSE (lines) are in agreement with MC 
simulations (symbols) for all cases; this verifies the reliability of the approximation scheme. 



V. STOCHASTIC RESONANCE 



Next, we study SR [39H48j in our model. SR is an intriguing phenomenon, and it plays 
an important role in systems accompanied by noise; hence, it has been studied extensively 
in various configurations. In particular, biological applications of SR have attracted con 



siderable attention, and they have been confirmed experimentally and theoretically 



49 



52| 



because biological mechanisms occur in noisy environments. We calculate the spectral am- 
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FIG. 4: (Color online) MFPT t s as a function of (a) the relaxation rate 7 and (b) the squared 
variation coefficient p. The lines and symbols denote results of CSE and MC, respectively, (a) 
D x = 1 and Q = 0.08 with p = 0.01 (solid line and circles), 0.1 (dotted line and squares), 1 (dot- 
dashed line and triangles), and 100 (dot-dot-dashed line and crosses), (b) D x = 1 and Q = 0.08 
with 7 = 0.3 (solid line and circles), 1 (dotted line and squares), 10 (dot-dashed line and triangles), 
and 100 (dot-dot-dashed line and crosses). The MC results are calculated as averages of 20000 
escape events. 

plification factor of SR with a periodic input under additive SIN. Specifically, we employ 
linear response approximation [53| to calculate the quantity. For a sufficiently small driving 
force, linear response approximation can be used to investigate SR. 

We assume that the system of interest is modulated by an external input eexp(— iflt), 
where e and Q are the input strength and the angular frequency, respectively. A Langevin 
equation is given by 

dx 

— = f{x) +eexp{-iQt) + sUt), (28) 
and Eq. (0), where f(x) = x - x 3 . The FPE of Eqs. (J2SJ and © is 

^P(x,s;t)=J? p P(x,s;t), (29) 

with 

d 

££ v = Jzf + JSfieexp(-ifit), Jgfi = -— , (30) 

ox 

where Jzfo is denned in Eq. (jSJ). We assume that e is sufficiently small for the system to be 
well approximated by the linear response. Let P as {x, s; t) be an asymptotic solution (t — > 00) 
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of Eq. (|29|) . According to the Floquet theorem, P as (x, s; t) is a periodic function having the 
same period as the input: 

P as (x,s;t) = P as (x,s;t + T), (31) 

where T is the period [T = 2tt/Q]. According to Eq. ( 13TT) and the linear response approxi- 
mation, we can expand P as (x, s; t) as 

P as (x 1 s;t) = P (x, s) + P\(x, s)eexp(— iQt). (32) 

From a normalization condition, Pi(x,s) must satisfy 

j dx j dsPi(x,s) = 0. (33) 

Substituting Eq. (1321 into Eq. (129]) and comparing the order of e, we obtain the following 
coupled equations: 

0(1) Sf P (x,s) = 0, (34) 
0(e) ^ Pi(x, s) + ^Poix, s) = -iQP^x, s). (35) 

Eq. fl34|) is identical to the equation for stationary distributions [Eq. ffl2l ]. Following the 
procedure for stationary distributions (Sec. IIIII) . we expand Pi(x, s) in terms of the orthonor- 
mal complete set. Using the relation P as {.x, s; t) = P as (—x, s;t + T/2) in Eq. ( 1301) . Pi(x, s) 
admits the odd symmetry expansion: 

K L 

Pi(x, s) = ip (x)4>o(s) ^2 G k/ Lp 2 k+i(x)ipi(s), (36) 

k=0 £=0 

where G k / are coefficients. Note that Eq. (136]) automatically satisfies Eq. (I3"3"|) because of 
the orthonormality. Following the same procedures as those in Sees. IIIII and IIV} Eqs. ([3] 
and (135]) can be represented as the following linear algebraic equation in terms of G k / 

= Cv / 2(2fc + l)C M + G k;i ^2k + l- -^{2k + l) 2 - 7 £ + ifij 



l-^ + 2( 2 D x ^a 2 + 2D sy } 

-Gk + i,e^^V(2k + 2)(2k + 3) - G? fc _ 2>/ ^2 VP + l)2fc(2A: - l)(2fc - 2) 
+2G^ lif+2 C 2 ^D sV ^pTIpdl(lT2) + 2G k - 1>e - 2 t 2 D x D 8 y/2k(2k + W ~ l ) 
+4G k -i,e-iC 2 D x a^2D s k(2k + 1)£ + AG k ^ +1 C, 2 D x a^2D s k{^ + !)(£+ 1). (37) 
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Ck,e has already been calculated in Eq. (|T8|) for the stationary distributions. 
From Eq. f )32|) . the time-dependent asymptotic average of x is given by 



( x ( t ))a S = j ' dx J dsxP as (x,s;t), 

= (x) + (x) 1 eexp(-m), (38) 



with 



<z> = / d X J ds ,*(,, .), («>,=/*: / <fe .), 

where (x) = owing to the symmetry. The susceptibility x is defined as the proportional 
coefficient of the input signal, which is given by \ = (x) 1 . There are several approaches to 
calculating the susceptibility e.g., the fluctuation-dissipation relation 36] or the moment 



method 



54 



55( | . Using the orthonormal and complete relations, the susceptibility is 



x = w (39) 

Let us consider a cosinusoidal input ecos{VLt). (x(t)) as for this case is 

( x (t)) as = (x) + \x\ecos(Qt + 6), 9 = - arctan (j^) , (40) 

where 9 is the phase. We evaluate the spectral amplification as \x\ 2 = |Go,o| 2 /(2C 2 )- 

We perform MC simulations to verify the reliability of the linear response approximation. 
For MC simulations, a method in Ref. [56 1 was employed. The averages of 2000 trajectories 
were calculated and the susceptibility was estimated by their variance [Eq. (|40p ] (the method 
of moments estimation). Fig. |5] shows (x(t)) as calculated by Eqs. ( 1391) and ()40l) (solid line) 
and MC simulations (circles). We observe excellent agreement between them, which verifies 
the reliability of the linear response approximation. 

Fig.|6]shows the dependence of the spectral amplification factor \x\ 2 on Q, p, and 7, where 
theoretical results obtained using CSE are denoted by lines and MC results are denoted 
by symbols. The MC results were in good agreement with those of CSE, verifying their 
reliability. 

Specifically, Fig. Eta) shows the \x\ 2 dependence on Q(a) = D X (D S + a 2 ) (a is varied 
while keeping D x and D s constant) with four 7 values: 7 = 0.3 (solid line and circles), 
7 = 1 (dotted line and squares), 7 = 10 (dot-dashed line and triangles), and 7 = 100 
(dot-dot-dashed line and crosses) with D x = 1, D s = 0.1, and Q = 0.1. Here, \x\ 2 achieves 
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FIG. 5: (Color online) (x(t)) as of CSE results [Eqs. fl39|) and (@DD] (solid line) and MC simulations 
as the average of 2000 trajectories (circles). The parameters are D x = 1, D s = 0.1, 7 = 1, a = 0.3, 
and n = 0.1, and e = 0.1 for MC. 

a maximum around Q(a) = 0.2, and the maximum is larger for larger 7. SIN approaches 
white noise for 7 — > 00, indicating that the strength of the SR effect is maximized under 
white noise. On the other hand, \x\ 2 in the range Q(a) > 0.4 has a different tendency, i.e., 
\x\ 2 is larger for smaller 7. Although the peaks of \x\ 2 at Q{a) ~ 0.2 are smaller for smaller 
7, SIN can induce better performance when the noise intensity exceeds Q(a) ~ 0.4. 

Next, we calculate \x\ 2 as a function of p with four 7 values: 7 = 0.3 (solid line and 
circles), 7 = 1 (dotted line and squares), 7 = 10 (dot-dashed line and triangles) and 7 = 100 
(dot-dot-dashed line and crosses). We vary p while keeping the effective noise intensity 
Q constant. Because the spectral amplification factor \x\ 2 is maximum as a function of 
the effective noise intensity in SR, its strength can be measured by the magnitude of the 
maxima. The maxima in Fig. EJ^a) are located around Q ~ 0.2; hence, we fixed Q = 0.2 
and investigated \x\ 2 dependence on p in Fig. fflb) (D x and Q are the same as as those in 
Fig. Eta)). Accordingly, \x\ 2 of Fig. Mb) can be identified as the strength of the SR effect as 
a function of p. Because SIN reduces to white noise as 7 — > 00, \x\ 2 as a function of p does 
not change for 7 = 100. On the other hand, \x\ 2 is more strongly affected by p for smaller 
7. SIN also reduces to white noise as p — > 0, and \x\ 2 increases as p — > in all cases. We 
observe non-monotonic behavior of \x\ 2 as a function of p, i.e., the strength of the SR effect 
is minimized around p ~ 1. Remarkably, the effect of the input signal is minimized around 



16 



(a) 



10 



Ixl 2 



r- 

- ■/ V. 

■ i ^ 

■ \ 


D x = 1 ■ 
D s = 0.1 - 

n = o.i - 







— theory 

7 = 0.3 

7 = 1 

7 = 10 

7 = 100 



- MC — 

o 7 = 0.3 

□ 7 = 1 

A 7 = 10 

x 7 = 100 





D x = 1 
Q =0.2 
Q = 0.1 



7 io 1 



FIG. 6: (Color online) (a) Spectral amplification factor \x\ 2 as a function of (a) the effective noise 
intensity Q(a) (vary a while keeping D s constant), (b) the square variation coefficient p, and 
(c) the relaxation rate 7. The lines and symbols denote CSE and MC results, respectively, (a) 
D x = 1, D s = 0.1 and ft, = 0.1 with 7 = 0.3 (solid line and circles), 1 (dotted line and squares), 10 
(dot-dashed line and triangles), and 100 (dot-dot-dashed line and crosses), (b) D x = 1, Q = 0.2 
and = 0.1 with 7 = 0.3 (solid line and circles), 1 (dotted line and squares), 10 (dot-dashed line 
and triangles), and 100 (dot-dot-dashed line and crosses), (c) D x = 1, Q = 0.2, and £1 = 0.1 with 
p = 0.01 (solid line and circles), 0.1 (dotted line and squares), 1 (dot-dashed line and triangles), and 
100 (dot-dot-dashed line and crosses). MC results with e = 0.1 are calculated using the variance 
of the average of 2000 trajectories (the method of moments estimation). 
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FIG. 7: (Color online) Illustrative description of the strength of the RA and SR effects as a function 
of the squared variation coefficient p. We observe a trade-off relation between the strength of RA 
and SR. 

p — 1, even though the strength of the noise intensity fluctuation is monotonic as a function 
of p. 

Fig. Etc) shows \x\ 2 as a function of 7 for four p values: p = 0.01 (solid line and circles), 
p = 0.1 (dotted-line and squares), p = 1 (dot-dashed line and triangles), and p = 100 (dot- 
dot-dashed line and crosses) with D x — 1, Q — 0.2, and Q = 0.1. In all cases, \x\ 2 increases 
as a function of 7; therefore, the SR effect achieves a maximum under white noise. 

VI. DISCUSSION 

In Sees. [IV] and |V] we have shown that the strength of the RA and SR effects exhibits 
non-monotonic behavior as a function of the squared variation coefficient p. Furthermore, 
the strength of RA and SR effects is enhanced in different p regions. The strength of the RA 
effect is maximum around p > 1, whereas that of the SR effect is stronger for p < 10~ 2 . On 
the other hand, the strength of the SR and RA effects is very weak in regions of p ~ 1 and 
p < 10~ 2 , respectively. These results shows that strength of these two effects has a trade- 
off relation in terms of p. An illustrative description of the trade-off relation between the 
strength of RA and SR effects is shown in Fig. [7j where the solid and dotted lines represent 
the strength of the RA and SR effects, respectively, as a function of p. 
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Langevin equations have been extensively applied to stochastic biochemical reactions 
such as gene expression [57| and neuronal response. In a zeroth-order approximation, these 



biological mechanisms can be modeled using a bistable potential [58] . Biological mechanisms 
are subject to many fluctuations having different time-scales. It has been reported theoreti- 
cally and experimentally that RA and SR are expected to play important roles in biological 
mechanisms. RA can minimize the delays in signal detection, which improves the response 
to signals. On the other hand, SR is responsible for accurate signal detection in noisy envi- 
ronments. These two factors are important in signal transmission, and our results indicate 
that their importance can be tuned with p. The results presented above may provide us 
with a new insight into the analyses of stochastic aspects of biological mechanisms. 

VII. CONCLUDING REMARKS 

In the present paper, we employed CSE to calculate stationary distributions, MFPT, and 



the spectral amplification factor. In our previous study [18(|, we used adiabatic elimination to 
derive a time evolution equation. CSE is advantageous in that the ranges of the relaxation- 
rate 7 and the noise intensity are not limited, as opposed to the adiabatic elimination-based 
method, which is valid for 7 > 1. In addition, CSE enables us to calculate quantities 
such as MFPT and the spectral amplification factor. On the other hand, using adiabatic 
elimination, we can calculate stationary distributions in the closed form, and it can be used 
for general non-linear drift terms. In contrast, CSE can only handle polynomial drift terms, 
for which stationary distributions are obtained by a numerical method. Both approaches 
are complementary. From the MFPT calculation, we identified the RA phenomenon as a 
function of 7. We also showed that the strength of the RA effect is highly dependent on the 
squared variation coefficient p, and that the strength of the SR effect as a function of p is 
minimum around p ~ 1. These results indicate that p, the ratio between the variance and 
mean of the noise intensity modulating process [Eq. f JTTj) ]. has a crucial impact on the RA 
and SR effects. 

Because CSE can be used for polynomial drift terms with arbitrary magnitudes of relax- 
ation rate and noise intensity, the analysis described in this paper can be applied to various 
real-world phenomena. Furthermore, we focused on periodic SR, in which the system of 
interest is modulated by a periodic input. With regard to biological cases, the investigation 
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of aperiodic SR [59|, |60[ is important. We plan to investigate this subject in the future. 
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Appendix A: Correlation function 

Here, we calculate the correlation function of SIN. By definition, the correlation function 
is given by 

(s(t)Ut)s(t')Ut')) = J dsds'd^dC [ss'U' x P(s^ x ;t\s',C^')P(s',C;t')]. (Al) 
Since s(t) and £ x (t) are independent, Eq. f lAlj) becomes 



(s(t)Ut)s(t')Ut')} = J dsds'd^ [ss'e^(s;<|s';0^(a , ;0^;<l^;< , )^;0], 

= (s(t)s(t')) (Ut)Ut')) , (A2) 
where the correlation function of s(t) is calculated as 

(s(t)s(t')} = D s exp {-j\t -t'\) + a 2 . (A3) 
From Eqs. flA2h and (IA30 . we obtain 



[s{t)Ut)s(t')Ut')) = 2D x {D s ex V (- 1 \t-t'\) + a 2 }5(t-t'), 

= 2D x (D s + a 2 )5(t-t'), (A4) 
= 2Q5(t-t'), (A5) 



where Q is the effective intensity defined by Eq. (J9J). From Eq. (IA5I) . the intensity of SIN 
is in agreement with the effective intensity Q, which is calculated via adiabatic elimination 
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